Adaptive Wavelet Collocation Method 
for Simulation of Time Dependent Maxwell's Equations 

Haojun Li*, Kirankumar R. Hiremath*'+, Andreas Rieder*, and Wolfgang Freude* 

^Department of Mathematics, Karlsruhe Institute of Technology, Germany 
+ Computational Nanooptics Group, Department of Numerical Analysis and Modelling 
Konrad-Zuse-Zentrum fiir Informationstechnik Berlin, Takustrasse 7, 14195 Berlin, Germany 
★Institute of High-Frequency and Quantum Electronics, 
Karlsruhe Institute of Technology, Germany 

Corresponding author: andreas.rieder®. kit.edu 
5 April, 2012 

Abstract 

This paper investigates an adaptive wavelet collocation time domain method for the numerical solu- 
tion of Maxwell's equations. In this method a computational grid is dynamically adapted at each time 
step by using the wavelet decomposition of the field at that time instant. In the regions where the fields 
are highly localized, the method assigns more grid points; and in the regions where the fields are sparse, 
there will be less grid points. On the adapted grid, update schemes with high spatial order and explicit 
time stepping are formulated. The method has high compression rate, which substantially reduces the 
computational cost allowing efficient use of computational resources. This adaptive wavelet collocation 
method is especially suitable for simulation of guided-wave optical devices. 

keyword: Maxwell's equations, time domain methods, wavelets, wavelet collocation method, adaptivity 

1 Introduction 

The numerical solution of Maxwell's equations is an active area of computational research. Typically, 
Maxwell's equations are solved either in the frequency domain or in the time domain, where each of these 
approaches has its own relative merits. We are specifically interested in efficient algorithms for light propa- 
gation problems in guided wave photonic applications [1], and work in the time domain. The most popular 
class of methods in this area is the finite difference time domain (FDTD) method [2]. Due to the structured 
grid requirement of these methods, they become cumbersome while dealing with optical devices having 
curved interfaces and different length scales. To overcome these difficulties, a discontinuous Galerkin time 
domain (DGTD) method has been investigated [3]. For a time dependent wave propagation problem, all 
these methods use a fixed grid/mesh for discretization. In general, such a grid can under-sample the tem- 
poral dynamics, or over-sample the field propagation causing high computational costs. If the spatial grid 
adapts itself according to the temporal evolution of the field, then the computational resources will be used 
much more efficiently. 

We propose an adaptive-grid method which represents propagating fields at each time step by a com- 
pressed wavelet decomposition, and which automatically adapts the computational mesh to the changing 
shape of the signal. In the initial studies of the wavelet formulation, the interpolating scaling functions 
were used for frequency domain waveguide analysis [4]. To the best of our knowledge, the suitability of 
the wavelet decompositions for time dependent Maxwell problems has not yet been investigated. Vasilyev 
and his co-authors developed the adaptive wavelet collocation time domain (AWC-TD) method as a general 
scheme to solve evolution equations, and they successfully verified the scheme's effectiveness in the area 
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of computational fluid dynamics [5, 6]. Based on these studies, we present in this work a proof-of-concept 
for an AWC-TD for the time dependent Maxwell's equations. 

The paper is organized as follows. In Sec. 2, we provide a brief account on Maxwell's equations and 
some of the related concepts for their numerical solutions. We start Sec. 3 with an introduction to (interpo- 
lating) wavelets, and how they can be used to discretize partial differential equations. Also in this section 
we explain the structure of AWC-TD method in the context of Maxwell's equations. Sec. 4 gives algorith- 
mic details of the method. Numerical results of the AWC-TD method are given in Sec. 5 which contains 
our numerical experiments of propagating a 2D Gaussian peak in homogeneous environment. Finally we 
close the paper with concluding remarks in Sec. 6. 



2 Time domain Maxwell's equations 

Propagation of optical waves in a linear, non-magnetic dielectric medium with no charges and currents is 
governed by the following time dependent Maxwell's equations 

-^^(f,t) = Vx^(f,i), ^^V{f,t) = V xn{r,t), V-P(f,t) =0, and V-i3(f,t) = 0, (1) 

where the electric field £ and the electric flux density V, as well as the magnetic field H and the magnetic 
flux density B, are related by the constitutive relations 

P(r,t) = eoer{r)£{f,t) and B{r,t) = noH{r,t). 

Here eq is the free space permittivity, is the relative permittivity and //q is free space permeability. 

For illustration purpose, we restrict ourselves to a 2D setting where the fields and the material properties 
are assumed to be invariant in the y-direction, i. e. r = {x, z) and the partial derivatives of all fields with 
respect to y vanish identically. We suppress the explicit function dependence on r and t. Then Maxwell's 
equations (1) decouple into a pair of independent sets of equations, 

d£^^ i_dHy d£^^j_&Hy my ^ r d£^ _ d£A 

dt eo£r dz ' dt eoSr dx ^ dt /iq \ dz dx J ^ 
identified as transverse electric (TE)y setting, and 

&H,^j_d£y &H^^_]^d£y d£y^j_rcm,_(mA 

dt fj,Q dz ' dt fio dx ' dt eo^r \ dx J ' 

identified as transverse magnetic (TM)j^ setting. Here £x, f^, • • • etc. denote the respective field compo- 
nents. 

Originally, Maxwell's equations are formulated for a whole space. For numerical computations we need 
to restrict them to a bounded computational domain Q as shown in Fig. 1 . This is done with a transparent 
boundary condition, which is reaUzed in our case with perfectly matched layer (PML) [7, 8]. The principle 
of PML is that (outgoing) waves scattered from the scatterer Jig pass through the interface between O and 
PML without reflections, and attenuate significantly inside the PML. The waves virtually vanish before 
reaching the outermost boundary of the PML, where the perfectly electric boundary (PEB) condition is 
employed. Implementation details about the PML technique specific for the method discussed in this paper 
can be found in Ref. [9]. For the sake of clarity, we work with the general formulation given by Eq. (2)-(3). 

As in the case of the standard FDTD method [2], in our approach we use the central difference scheme 
for the time derivatives in Eq. (2)-(3), but we will construct a different discretization scheme of the spatial 
derivatives. This is done with interpolating scaUng functions and lifted interpolating wavelets (explained 
in Sec. 3). The induced multiresolution approximation [10, 11] enables us to decompose fields into various 
resolution levels, and thus allows to discard unimportant features. As a result, we will obtain a variant of the 
FDTD method, which is constructed with respect to a locally refined grid. In the next section we describe 
this numerical scheme in detail. 
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Figure 1 : Typical simulation setting with a computational domain J7 surrounded by the perfectly matched 
layer. Here just for the sake of illustration, we show the scatterer fig completely enclosed inside Other 
configurations like incoming-outgoing waveguides are also possible [9]. 



3 Adaptive wavelet collocation method 

The adaptive wavelet collocation (AWC) method was proposed by Vasilyev and co-authors in a series of 
papers [12, 13, 5, 6] as a general scheme to solve evolution equations. In the present section, we tailor the 
AWC method to tackle Eq.(2)-(3). In contrast to the originally formulated AWC method, we do not need to 
utihze second generation wavelets, which have been mainly invented to implement boundary constraints, 
and to find wavelet decompositions on irregular domains. Since we use the PML method, we can identify 
field values outside the PML region with zero, and therefore we are not forced to adapt our wavelets to 
the boundary restrictions. Hence, we consider only the first generation wavelets, which are generated by 
the shifts and the dilations of a single function. Now we outline the essential steps for computing spatial 
derivatives of functions in wavelet representations. 

3.1 Preliminaries 

A starting point of the AWC method is a wavelet decomposition of a function / G L^(M): 

+ 00 



where jo G Z, is the scaling function and ip is the wavelet function [14, 15]. For all j,n G Z, by (f)j^n 
and ipj,n we abbreviate the dilated and translated versions of 4> and ip, i.e. ^j,n(-) = 2^/'^4>{2^ ■ -n). 



The first (single) sum in (4) represents rough or low frequency information of /, while the second 
(double) sum contains the detail information at various resolution levels starting from the level jo to -|-oo. 
The absolute magnitude of the coefficients ajg^A; and jij^m measure the contributions of (pj^^k and i^j^m 
to /. By discarding terms in the double sum for which the wavelet coefficients /3j.m are absolutely less 
than a given threshold, one can efficiently compress the representation of /. This wavelet decomposition 
compression principle is exploited in the AWC method to enhance the computational efficiency. 

There are various families of the scaling functions (j) and wavelet functions ^ allowing representations 
like (4). As in [5, 6], we work with the interpolating scaling functions [16] and the corresponding lifted 
interpolating wavelets [17, 18]. Due to their interpolation property, we have 



can be seen as a variant of the well known finite difference method. We exploit this interpolating property 
in Sec. 3.2 and Sec. 3.5. 




(4) 
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In particular, we use the interpolating scaling function (ISF) family developed by Deslauriers and 
Dubuc [19, 16]. They constructed the interpolating functions by the iterative interpolation method, which 
does not require the concept of wavelets. Later Sweldens [17, 18] constructed the corresponding wavelet 
by lifting the Donoho wavelet [20]. We use DDn to denote ISF of order N, and L'^^ to denote the lifted 
interpolating wavelet of order A'^. Here the order N means that any polynomial p of degree k < 2N — 1 
can be expressed as 

p(-) = ^ CmDDNi- - m) 

m 

with suitable coefficients {c^}. The order N is half the number of the vanishing moments of the Ufted 
interpolating wavelet, i.e.. 



/ 



x''Dl^{x)dx = 0, k = 0,l,...,2N -1. 



Further details can be found in [17, 18, 9]. We normally choose same orders for the ISF and the lifted 
interpolating wavelet, i.e., N = N. It is easy to see that DDn and 01^;^ have compact supports, which 
increase with the order N. 

For the TMj^ setting in Eq. (3), the electric and magnetic fields depend on the spatial variables {x, z). 
As usual, see, e.g., [11, 15], we represent 2D fields by expansions of 2D scaling functions and wavelets 
which are defined by 



(f)Nix,z) := DDn{x)DDn(z), 

' DIn{x)DDn{z) 



DDn{x)DIn{z) 
[ DIn{x)DIn{z) 



= 3, 



and use the following abbreviations 

{(t)N)j,m,n{x,z) := {DDN)j^^{x){DDN)j^niz), 



{ (D/;v),,„(x)(DDiv),+i,2n(^) 



iDDN),+,,2mi^)iDlM)j,^{z) 
[ iDlN)j,mi^mN)j,Jz) 



1^ = 3. 



Let jmin and jmax (with jmin < jmax) be the coarsest and the finest spatial resolution levels. Let us consider 
/ G L^(M^) with exact resolution level jmax, that is, 

/ ~ ^ ^ '^jmax,m,n(0jv)jmax,m,n- (5) 

Then the wavelet representation of / with coarsest resolution level jmin is given by 

3 Jmax 1 

/ = Q'imin,m,n(0Ar)jinin,m,n + Pi,m,n{'^N)j,m,n (6) 

where the scaling coefficients {oij^^^,rn,n} and the wavelet coefficients {/SJ^ „} can be calculated from the 
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level jmax scaling coefficients {cej^^^rn,n} by the normalized 2D forward wavelet transform (FWT): 

dj,m,n = 2 (cj+l,2m+l,2n - ^ 2s_;Cj+i,2m+2i,2n) , (7a) 

I 

d'j,m,n = 2 (ci+l,2m,2n+l - ^ 2s_iCj+i,2m,2n+2i) , (7b) 

I 

^j,m,n = ^ {cj+i^2m+l,2n+l " ^ 2s_iCj+i_2m+2i,2n+l " ^ 2s_;/Cj+i_2m+l,2n+2F 

+ XI Xl*^^^-')*^^^-^')'^-?+l'2"»+2',2n+2i') , (7C) 
« I' 

Cj,m,n = Cj+l,2m,2n + ^ S-i4,"i+',n + XI ^-l'^lm,n+l' + X XI ^-l^-l"^%m+l,n+l' ^ (7d) 



with the following normalization conventions 

The coefficients 2s; and are Lagrangian interpolation weights. For example, when N = 2, these weights 
are 

s_2 = -l/16, s_i = 9/16, So = 9/16, si = -1/16, 2s_i = -1/16, 

and 

2so = 9/16, 2si = 9/16, 2s2 = -1/16. 

Readers may consult [16, 21] and [17, Theorem 12] for an explanation of how and why Lagrangian weights 
enter the iterative interpolation process. 

We also can compute back from the wavelet representation (6) to the scaUng function representation (5) 
by the inverse wavelet transform (IWT): 



Cj+l,2m,2n = Cj,m,n - X ^-l^],m+l,n + X ^-l"^],m,n+l' + X X *-'^-''4,m+i,n+/'' (^a) 
I V II' 

Cj+l,2m+l,2n = 2c^j,m,„ + X 2S-/Cj+i,2m+2i,2n, (8b) 

I 

Cj+l,2m,2n+l = '^d^j,m,n + X 2S-iCj+l,2m,2n+2i, (8c) 

I 

Cj+l,2m+l,2n+l = ^d^j,m,n + X 2^-''^J+l,2m+2;,2n+l + X 2^-'"^j+l,2m+l,2n+2i' 

/ V 

- y^^y^,{'^S-l){2S-li)Cj+i^2m+2l,2n+2l'- (8d) 
I I' 

3.2 Adaptive grid refinement wavelet compression 

We thin out the triple sum in (6) by discarding small wavelet coefficients, which corresponds to small scale 
details. For a given threshold C > 0, let 

3 Jmax 1 

/C X ":'min,m,n(0jv)jmm,m,n + X X X '^C ^^j,m,n)i''pN)j,m,n, 

m,n u=l j=jmin m,n 



where the threshold function : M ^ M is defined by 



X 





for u e {1, 2} and |x| > 2-^-V2^^ 
for = 3 and \x\ > 2~^(, 
otherwise. 
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Note that we have defined the uniform threshold C in terms of the normahzed wavelet coefficients d^^ ^ 
defined in Eq. (7) i.e., if „| < C then dj ^^^^ = in in Eq. (6). Then the compression error is 
proportional to ( [5]: 

||/-/clloo<cc. 

Our basis functions in (6), which are translates and dilates of </)jv and ij^'^, are interpolating at the corre- 
sponding grid points. Let 

m n 
Xj,m ■= ^ and Zj^n ■= ^ for m,neZ, 

then we have the following one-to-one correspondence between the basis functions and the grid points: 

j,m,n ^ (^j+l,2m5 Zj+l,2n+l)i iM^ N)j,m,n ^ (^j+l,2m+l; Zj+l,2n+l)- 

Here this correspondence means the validity of the interpolation property. For instance, we have that 

i.4^N)j,m,n{Xj,m'^ Zj,n') ~ ^m,m'^n,n' • 

With this explanations, we justified the synonymous usage of compression of the wavelet representation 
and compression/adaption of the grid points. 

3.3 Adjacent zone 

With the above described wavelet compression, the grid gets suitably sampled only for the current state 
of the fields. For a meaningful (i.e. physical) field evolution in the next time-step, the grid need to be 
supplemented by additional grid points, on which the fields may become significant in the next time step. 
This allows the grid to capture correctly the propagation of a wave. To this end Vasilyev [5, 6] has introduced 
a concept of an adjacent zone. 

To each point P = {xj^m, -2^j,n) in the current grid, we attach an adjacent zone which is defined as the 
set of points (a;j',m') -2j',n') which satisfy 

\f - j\ < L, \2^'-^m - m'\ < M, \2^'-^n - n'\ < M, 

where L is the width of the adjacent levels and M is the width of the physical space. As in [5], we verified 
that L = M = 1 is a computationally sufficient choice. Then the adjacent zone for a point P can be 
depicted as in Fig. 2. 

® neighbor points of the same level 

® 

• neighbor points of one finer level 

© 

Figure 2: Description of the adjacent zone of a grid point P. 

Note that the concept of adjacent zone is reasonable only for continuously propagating waves, as in case 
of our guided-wave applications, where in each time step the propagating waves do not travel far from the 
current position due to their finite propagation speed. 

3.4 Reconstruction check 

In this work we use the wavelet decompositions of the fields only to determine the adaptive grid. We do 
not propagate fields in their wavelet representations (cf. the statement in the first paragraph of Sec. 4). Thus 
at each time step, after adapting the grid using the FWT, and adding the adjacent zone, we need to restore 
the fields in the physical space by performing the inverse wavelet transformation (IWT). To this end, we 



• OP • 

• • • 
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may need to augment the adaptive grid with additional neighboring points (e.g. see Fig 3). This process 
of adding neighboring points needed to calculate the wavelet coefficients in the next time step is called 
reconstruction check. Fig 3 shows various possible scenarios, and the corresponding minimal set of the grid 
points required for calculation of the wavelet coefficients. The values of the wavelet coefficients at these 
newly added points are set to zero. 
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Xj,m-1 Xj,m Xj,m+1 Xj,m+2 

(b) x: The point correspond- 
ing to d|,^,„; •: The neigh- 
boring points needed to calculate 



yj-"+2— ^ — i ::[) i — i— 

yj.m+1 - 



yj,m 



yi,m-l 



t t ' t t 

Xj,m-1 Xi,m Xi,m+1 Xj,in+2 



(c) X : The point corresponding 
to rff.m,™; °: The neigh- 
boring points needed to calculate 



Figure 3: Descriptions of the neighboring points needed to calculate the wavelet coefficients d 
orders N = N = 2. 



with the 



The efficiency of the wavelet transform depends on the number of the finest grid points only at the 
beginiung; however, after the first compression, it depends solely on the cardinahty (= number of grid 
points) of the adaptive grid. 

3.5 Calculation of the spatial derivatives on the adaptive grid 

After the adjacent zone correction and the reconstruction check, we are in a position to calculate the deriva- 
tive of at a grid point in the adaptive grid. For this we need to know the density level of this point, which 
is defined as the maximum of the x-level and the z-level of that point. 

We illustrate this concept explicitly only for the x-level, the 2;-level can be determined analogously. For 
a point Q = (xq, ^o) in the adaptive grid Q, let Q' = {xi,zo) G ^ be the nearest point to Q. Then the 
x-level Levelx of Q relative to Q is 



Levelx := jmax - log2(dist((5, Q')/Ax) 



(9) 



where Ax is the smallest computational mesh size along the x axis, and dist((5, Q') = |xi — xo|. For 
dist(Q, Q') = Ax, the level Levelx of Q attains its maximum jmax- For dist((5, Q') = 2Ax, we have 
Levelx = jmax — 1> etc. See Fig. 4 for an example of describing the density level of a grid point. 



N 
<l 



2Ax 



Figure 4: Description of the density level of a point Q in an adaptive grid: the x-level of Q is jmax — 1 and 
the 2;-level of Q is jmax, thus, the density level of Q is jmax- 



Now we continue to discuss the derivative calculations. Suppose jo to be the density level of Q in ^. 
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Then, we can represent by a finite sum Pj^f locally in some neighborhood of Q. 

i<l>N)jo,m,n{x,z), (x, z) G (10) 

m,n 

We differentiate Pjo/ with respect to x to approximate the x-derivative of / at Q. If any points in the sum 
(10) are not present in Q, then we interpolate the values at these points by the IWT using the values of the 
coarser levels. From the interpolation property of {4>N)jo,m,n we know that 

ajo,m,n = 2-^°{PjJ) (^^, ^) , for m,neZ. 

Thus, we have 

{PjJ){x, z) = J2(PjJ) ^)dD^{2^°x - m)DDj,{2^^z - n), (x, z) G Oq. (11) 

m,n 

Differentiate both sides of (11) with respect to x gives 

m,n 

The derivatives of DD'^j can be calculated exactly at the integers using the difference filters shown in Table 1 
(see Ref. [16] for details of the derivation). 



i 


iV = 2 


= 3 


iV = 4 


1 


2/3 


272/365 


39296/49553 


2 


-1/12 


-53/365 


-76113/396424 


3 




16/1095 


1664/49553 


4 




1/2920 


-2645/1189272 


5 






-128/743295 
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1/1189272 



Table 1: Difference filters {DD'j^{i)}i^i with consistency order 2N. Note that DD'^{i) = -DD'^{-i). 



Since the density level of Q is jo> there exist m', n' ^"L such that Q = ^) and it is easy to see 

that 

^^fe'^j -2^(Pio/)(^,^J 5^ ^^iv(n -n) 

m,n 

= 2-^(P,J)g,|^)l)I)^K-m). (13) 



Similarly, 



This finishes the general discussion about the adaptive wavelet collocation method; in the next section, we 
apply it to Maxwell's equations. 



4 AWC-TD method for Maxwell's equations 

In this section we formulate the update scheme for Maxwell's equations, and then elaborate on algorithmic 
issues related with the AWC-TD method. In the present formulation we represent the electric and magnetic 
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fields in the physical space, and not in the wavelet space. To unleash the full power of adaptivity, however, 
the field representation and the update in wavelet space are advantageous. 

We illustrate the method for the transverse magnetic (TM)^ setting given by (3). Similar procedure 
can also be formulated for TE^ setting in (2). Unlike the standard FDTD method, here the electric field 
and the magnetic field components are evaluated on same spatial grid, and their spatial derivatives are 
approximated at the same grid point. But the electric field components are sampled at integer time-steps, 
whereas the magnetic field components are sampled at half-integer time-steps. 



4.1 Update scheme for the spatial derivative 

For a point Q in the adapted grid Q, let T-Lxtq^^"^ , ^zIq^^^^ ^'^'^ ^y\^ denote the discretized value of 
Tiz and £y at the point Q, and at a time (A; + 1/2) At for the magnetic field components and at a time A; At 
for the electric field component where At > is the time step size (Note that, the electric field components 
are sampled at integer time-steps, whereas the magnetic field components are sampled at half-integer time- 
steps.). Assume j{Q) to be the density level of Q relative to Q. Then we can represent the point Q as 
{xj{Q),m',Zj{Q),n') for somc m',n' ^ Z. 

Let L be the length of the computational domain Q,. We rescale the wavelet decomposition (11) with 
the factor L. Then using the central difference scheme for the time derivatives and using (13)-(14) for the 
spatial derivatives, we get the following difference equations 

^ n 

h^i k-i At2''('5).^ 
nSn^=nAn^ + — y^£v\f-r r ^DDUm'-m), (15b) 

^ m 

""^ '"'^ So Erlg L \^ '(^j(Q),m'>^j{Q),n) ' 

-Y'^ztt^^ , .L>L>^(m'-m)V (15c) 

m 

The first time step {k = 0) is an expUcit Euler step with step size At/2 using initial conditions for the fields 
at the time t = 0. If not expUcitly mentioned, otherwise the fields are set zero at the beginning for all our 
numerical experiments in Sec. 5. The update equations for the PML assisted Maxwell's equations can be 
found in Ref. [9]. 

From the form of these update equations, it is clear that the AWC-TD method can be thought as an 
variant of high order FDTD method. The AWC-TD method is defined with respect to a locally adapted 

mesh, and unlike the FDTD method, it does not require a static (fixed), structured mesh. This will lead to 
efficient use of the computational resources. In the next section, we elaborate on algorithmic aspects of the 
method. 



4.2 Update scheme for the time derivative 

Several choices are available for time stepping. As in case of the standard FDTD method, we use in (15) 
the central difference scheme for the discretization of the time derivatives. For this explicit scheme, the 
smallest spatial step-size restricts the maximal time-step according to the Courant-Friedrichs-Lewy (CFL) 
stability condition. Using a uniform spatial mesh in the update equations (15) with a mesh size A in both 
coordinate directions the CFL condition reads 



A 



^t<— p-j , (16) 



see [9, Sec. 3.5] and [22], where c is the speed of light in vacuum and {DD'j^{l)} is the known derivative 
filter of the ISF as in Table 1 . Due to the local adaptive grid strategy of the AWC-TD method, we cannot 
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define a global stability criteria as above. But choosing A to be the smallest step size in the adaptive grid, 
we get a conservative bound for At via (16) for the AWC-TD method. In the simulation tests (in this paper, 
and in [9]) we did not experience any stability related issues with this modus operandi. 

4.3 Implementation aspects 

4.3.1 Grid management 

In AWC method the computational grid is changed with the state (spatial localization) of the propagating 
field. Thus the grid management is one of the important steps in the implementation of this method. This 
is done as following: We store the information of the adaptive grid into a 2D Boolean array called a grid 
mask or simply a mask, whose size is square the number of the finest grid points along one direction. We 
use 2D arrays of real numbers with the size of the grid mask to store the fields such as £y, Hx, and Tiz etc. 
Note that the computational effort for updating the fields at each time step is proportional to the cardinality 
(i.e. the number of entries in the mask with value 1) of the adaptive grid. 

If the value of an entry of a mask is true or 1, then the corresponding grid point is included in the 
adaptive grid; otherwise, it is not included in the grid. Thus by forcing the value of an entry of a mask 
to 1, we can include the corresponding point to the grid, or by forcing the entry to 0, we can exclude the 
corresponding point from the grid. 

4.3.2 Algoritlimic procedures 

Algorithm 1 outlines the main function awcm_main() of AWC-TD method for TMj^ setting. It mainly 
consists of two blocks of operations: The first block is initialization, and the second block is time stepping. 
In the time stepping block, at each time step the routines awcm_adaptive( ) and awcmMpdate( ) are called. 
The former routine optimally adapts the computational grid for the field updates at the next time step, 
whereas the latter routine calculates the spatial derivatives on the non-equidistant, adaptive grid, and updates 
the field values. 



Algorithm 1: awcm_main() for TM^^ settings 

# Initialization 
awcmJnitializeO 

# 

# Time stepping of £y, T-Lx and Tiz 
for t < T do 

# Adapt the grid for t + At according to £y, see Algorithm 2. 
awcm_adaptive() 

# 

# Update n*,'^^*^^, nl'^^^^^ and Sj+^K see Algorithm 3. 
awcm_update() 

# 

# Go to the next time step. 
_ t = t + At 



The initialization subroutine awcm_initialize() ensures that various required inputs for the AWC method 
are systematically prepared. It consists of checking the given initial data (i.e. for time step k = 0) Hx'"^, 
T-Lz^^ and £y^ at the finest resolution level jmax, the threshold (, the maximum and the minimum spatial 
resolution levels jmax and jmin respectively, and the number of time steps kmax- The time step At is chosen 
such that it satisfies the CFL condition given by (16). 

The adaptivity procedure in Algorithm 1 handled by a subroutine awcm_adaptive() is outlined in Algo- 
rithm 2. It is done by means of a 2D array £y with a mask MaskO. For later use, we store a copy of MaskO 
in pMaskO, since M askO will be modified by the subsequent subroutines. The duplicate pMaskO serves 
as a reference for finding those points which need to be interpolated before we can update the fields. We 
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perform the fast wavelet transform of Sy on M askO. Note that M askO is either fully 1 (as at the beginning) 
or a reconstruction check has been performed in the previous time step. In any case, FWTs on M askO are 
always possible. By the FWT applied to £y we obtain the scaling coefficients on the coarsest level jmin> 
and the wavelet coefficients on levels from jmin to jmax — 1- 

For each wavelet coefficient, we compare its absolute value with the given tolerance C- If it is less than 
C, we remove the corresponding point from MaskO. Next, we determine the adjacent zone for each point 
in MaskO, and then modify MaskO to include all points in these adjacent zones. Finally, a reconstruction 
check is applied to MaskO so that the FWT in the next time step is well defined. The latter two processes 
are done in the subroutine Maskext(Mas/i;0) as shown in Algorithm 2. 



Algorithm 2: awcm_adaptive() for TMj^ settings 

# Store MaskO into pM askO 

# pMaskO: The adaptive grid for ty at current time step. 
pMaskO = MaskO 

# 

# Fast wavelet transform of £y on MaskO with 

# £y is converted into coefficients of wavelet domain, M askO is thinned. 
FWT(Sy, MaskO, Q 

# 

# Add adjacent zone and perform a reconstruction check to MaskO. 
Maskext(MosfcO) 

# 

# Add points needed to calculate and -q^ on MaskO. 

# 1. Determine the density level of each point in MaskO. 
LevelO = Le\el(MaskO) 

# 2. Initialize Maskl with MaskO. 
Maskl = MaskO 

#3. Update Maskl. 
gMaskext(MasA;l, LevelO) 

# 

# Add points needed to calculate ^^j^ and on Maskl. 

# 1. Determine the density level of each point in Alaskl. 
Levell = Level(MasA:l) 

# 2. Initialize Mask2 with Maskl. 
Mask2 = Maskl 

#3. Update Mask2. 
gMaskext(MasA;2, Levell) 

# 

# Inverse wavelet transform of the values £y in the wavelet domain on Mask2. 

# £y is reconstructed from the values in the wavelet domain on M ask2. 
1WT(£:„, Mask2) 



After the above adaptation of the grid is done, we still need to make further reconstructions on this grid, 
so that it will allow com|)utation of the field derivatives required for the field update. For updating T-Lx and 
Tiz, we need -q^ and -g^ (see (3) or (15)). To calculate these spatial derivatives of the electric field, we 
interpolate values of £y at those neighbors of points in MaskO which are not already in MaskO. We store 
the information of MaskO into Maskl. Further, we add all points to Maskl needed in the calculations 
of spatial derivatives according to the density levels of the points in Maskl. These density levels are 
computed in subroutine Level(MasA;l) and stored in the 2D array LevelO. Again a reconstruction check of 
Maskl is required to enable IWTs. This is done by the subroutine gMaskext(MasA;l, LevelO). 

Then we need to follow the same procedure as above for updating £y using the spatial derivatives 
and ^2^. Again we add the neighboring points needed for calculations of the spatial derivatives of the 
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Algorithm 3: awcm_update() for TM^ settings 



# io update H^'- 

# Interpolate l-Li on points in Maskl which are not in pMaskO using inverse wavelet transform. 
interpolateCHa;, pMaskO, Maskl) 

# Calculate ^ on Maskl using Algorithm 4. 
dAz = diflz(£y, Maskl, Levell, dfilter, dz) 

Update n X on Maskl using dA^ as per formulation in (15a). 

# 

# To update Tiz- 

# Interpolate V.^ on points in Maskl which are not in pM askO using inverse wavelet transform. 
interpolateCW^, pMaskO, Maskl) 

d£ 

# Calculate -g^ on Maskl using Algorithm 4. 
dAx = diffx(£^j,, Maskl, Levell, dfilter, dx) 

Update Hz on Maskl using dAx as per formulation in (15b). 

# 

# To update £y: 

# Interpolate £y on points in MaskO which are not in pMaskO using inverse wavelet transform. 
interpolate(f^y, pMaskO, MaskO) 

# Calculate and ^ on MaskO, see the Algorithm 4. 

# diffx() is defined in Algorithm 4. diffz() is similarly defined. 
dAz = diffz(?^^, MaskO, LevelO, dfilter, dz) 

dAx = diffxCHz, MaskO, LevelO, dfilter, dx) 

Update £y on MaskO using dAz and dAx as per formulation in (15c). 



magnetic field. We copy Maskl to Mask2, and calculate the density level array Levell of Mask2. The 
necessary reconstruction check is then done by calling gMaskext(MasA;2,Lef e/1). The call of JWT{£y, 
MaskT) to reconstruct £y in the physical domain finishes the routine awcm_adaptive() in Algorithm 2. 

Next, we update the field values on the adaptive grid, which is described by Algorithms 3. Since the 
adaptive grid may change with time, we need to interpolate the field values at points in the adaptive grid of 
the current time step, which are not included in the adaptive grid of the previous time step. For example, 
consider the update of Hx about a grid point Q at a time {k + 1/2) At in (15a). Since Q is not necessarily 
in the adaptive grid of previous time {k — 1/2) At, the value T-LxIq ^^'^ in (15a) must be interpolated. Once 
this is done. Algorithm 4 calculates the spatial derivatives of each field components on the adaptive grid, 
and then the fields are updated. 

5 Numerical results: Gaussian pulse propagation 

In this section we demonstrate the applicability of the AWC-TD method. The method has been implemented 
in C++, and the computations have been performed on 32 GB RAM, Linux system with AMD Opteron 
processors. 

As an example, we consider propagation of a spatial Gaussian pulse in free space (e^ = 1)- We solve a 
system of TMj, equations within a square domain i7 = [— L/2, L/2] x [— L/2, L/2] in the XZ plane. We 
set the domain length L = 6.0 fim, the PML width d = L/4, and the initial spatial Gaussian excitation 
£y{x, z, 0) = exp(— (x^ + z'^)/{2a'^)) with the Gaussian pulse width a = 1/(4 \/2) /um. Implementation 
details about the PML can be found in Ref. [9]. 

Our minimum and maximum resolution levels are jmm = 3 and jmax = 9 inducing the smallest mesh 
size A = Ax = Az = L/2-''""^ = 11.71875 nm. The temporal error of the AWC-TD method is controlled 
by O(At^) if we do not consider the compression, which is the consistency order of the central difference 
discretization of the time derivatives. Accordingly, a reasonable choice for the threshold ^ is a value slightly 
larger than the discretization error. As the orders of the underlying interpolating scaling function/wavelet 
pair is = iV = 4, we set At = A/c/1.6, which is just below the maximal step size from the CFL 
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Algorithm 4: diffx(A, Mask, Level, dfilter, dx) 



Input : A = 2D array of field values 
Mask = grid mask, 

Level = x-level of each point in Mask, 
dfilter difference filters as given in Table 1 , 

dx = the smallest mesh size in x direction at the highest resolution level 
Return: a 2D array of ^ on Mask 
# Initialize a 2D array dA for the storage of 

dA = 

M = {{^j,m,yj,n) I m, n = 0, 1, . . . , 2-'}, where Xj^rn = -^,yj,n = for jmin <j< jmax- 

forall the Q = G /Cj„,,, do 

if Q G Mask then 

# Read the density level of Q from Level. 
j{Q) = -Z^eve/ [n] [m] 

Calculate dA at point Q using dfilter and values of A at neighbor points in the level j{Q) 
\_ as described in (13). 



condition (16). For this setting, a choice of wavelet threshold ( = 5.0 x 10 experimentally turned out 
to be sufficient concerning both adaptivity and accuracy. The Gaussian pulse, launched in the center of the 



t=20At, cp=0,88% t=200At, cp=1.41% t=380At, cp=2,02% t=560At, cp=2,40% 




t=740At, cp=2.73% t=920At, cp= 1 .76% t= 1 lOOAt, cp=0,80% t= 1280At, cp=0,61% 




Figure 5: Evolution of the initial excitation £y{x, z, 0) = exp(— (x^ + z'^)/{2a'^)) in the XZ plane with 
o" = 1/(4 \/2) fim, N = N = 4 and ( = 5.0 x 10^^. On top of each time frame, time and grid compression 
rate cp are given (cp is the ratio of the cardinality of the adaptive grid and the cardinality of the full grid 
with a uniform step size A (= the smallest mesh size) in the both coordinate directions). The adaptive 
grid systematically follows and resolves the wavefront. In regions where the field is small or not present 
only grid points of the coarsest level are assigned. For an animation movie, see the YouTube channel: 
www . youtube . com/user/Hao junLi#p/u/ 1 /2Yzp j f 7Xnp4. 

computational domain, spreads away from the center as time evolves. Fig. 5 illustrates how the adaptive grid 
systematically follows and resolves the wavefront. Since the electromagnetic field energy is spreading in all 
directions, the field's amplitude is decreasing (unlike as in ID, where during the propagation the amplitude 
stays at half of the initial value, see [9, Sec. 4.4.1]). The AWC method generates a detailed mesh only in the 
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regions where the field is localized, the mesh gets coarse in other parts of the computational domain. As 
seen in the snapshots for t = 200At or t = 920At, it is evident that depending on the extend of the field 
localization, the density of the grid points varies accordingly. 

A figure of merit for the performance of the AWC-TD method is the compression rate cp, which is 
defined as a ratio of the cardinality of the adaptive grid and the cardinality of the full grid with a uniform 
step size A (= the smallest mesh size) in the both coordinate directions. The percentage cp on the top of 
each time frame in Fig. 5 shows the grid compression rate. Since the extent of a spatial localization of a 
pulse depends on its frequency contents, the compression rate cp for the test case in Fig. 5 varies (also seen 
in Fig. 7). Nevertheless, for all time steps the number of grid points in the adapted grid is substantially less 
than that of in the full grid; but still the AWC method resolves the pulse very well with an optimal (with 
respect to the given threshold Q allocation of the grid points. 
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Figure 6: Relative Error in £y between the adaptive wavelet collocation method and the full grid wavelet 
method. 

The relative maximal error of £y field values over Q between the adaptive and the full grid methods as 
the time evolves is shown in Fig. 6. Despite of grid compression (which can be quite significant at some 
time instants, as seen in Fig. 5), the solution by the AWC method is quite close to that of by the full grid 
method. As mentioned earlier, as the pulse spreads in all the direction, the field becomes weak, and the real 
performance gain by the adaptivity effectively reduces. It is reflected in the apparent increase in the relative 
maximal error (with respect to the full grid method) in Fig. 6. Note that when the field has completely left 
the computational domain 0, roughly after 800 time steps, the error over Q is not defined meaningfully any 
more. 

Fig. 7 demonstrates that (the major part of) the computational effort of the AWC-TD method per time 
step is indeed proportional to the cardinality of the adapted grid at that time instant. To this end, we recorded 
the CPU time for every ten time steps (Fig. 7 top). For comparison, we also plotted the grid compression 
rate as a function of the time step (Fig. 7 bottom). Both functions progress in parallel, thus validating the 
above assertion about the numerical effort of the AWC-TD method. 

6 Conclusions 

In this paper we investigated an adaptive wavelet collocation time domain method for the numerical solution 
of Maxwell's equations. In this method a computational grid is dynamically adapted at each time step by 
using the wavelet decomposition of the field at that time instant. With additional amendments (e.g. adjacent 
zone corrections, reconstruction check, etc.) to the adapted grid, we formulated explicit time stepping 
update scheme for the field evolution, which is a variant of high order FDTD method, and is defined with 
respect to the locally adapted mesh. We illustrated that the AWC-TD method has high compression rate. 
Since (the major part of) the computational cost of the method per time step is proportional to the cardinality 
of the adapted grid at that time instant, it allows efficient use of computational resources. 

This method is especially suitable for simulation of guided-wave phenomena as in the case of integrated 
optics devices. Initial studies for simulation of integrated optics microring resonators can be found in [9]. 
In the present feasibility study we represented the electric and magnetic fields in the physical space, and 
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Figure 7: CPU time (top) and grid compression rate cp (bottom) as functions of the time step. Both functions 
progress in parallel which illustrates the fact that the numerical effort of the AWC-TD method for each time 
step is proportional to the number of points in the actual grid. 

not in the wavelet space. To unleash the full power of adaptivity, however, the field representation and the 
update in wavelet space are mandatory. 
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